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I. INTRODUCTION 

Since its first appearance in the Density Matrix 

Renormalization Group methodoa has witnessed great 
developments and it soon became one of the most widely 
applied numerical methods in one-dimensional solid state 
physics. Within a short period of time, the real space 
renormalization method had been further extended and 
the momentum space version of the method (fc-DMRG) 
was introduced by Xiangta in 1996. Unfortunately, test 
calculations on the Hubbard model indicated relatively 
poor performance compared to the real space version 
which hindered further application of the method for sev- 
eral years. 

Quite recently, DMRG was used to study mod- 
els of cyclic polyenes!] and models of polyacctylcncEl. 
S. R. White has successfully applied fc-DMRG in quan- 
tum chemistry to calculate the ground state energy of 
molecules represented in the framework of the usual Lin- 
ear Combination of Atomic Ophitals (LCAO) approxi- 
mation, using small basis sets&Q. His results seemed 
challenging and attracted-,|eonsiderable attention which 
stimulated various groupsocl to start to work on the new 
field. 

Among all the various models studied by DMRG dur- 
ing the past decade the accuracy of the algorithm has 
always been a problem which is still not satisfactorily 
solved. The recent application of DMRG in quantum 
chemistry gives further grounds for benchmark investi- 
gations of this question within the new framework. In 
all attempts so far, the accuracy of the method was ana- 
lyzed a posteriori by means of comparison with the corre- 
sponding full CI (FCJ) benchmark results. For instance, 
recently, Chan et aln reexamined the scaling behavior of 
the real error, developing an extrapolation approach as a 
function of the number of block states (M) . 

In this paper we show that in contrast to previous ap- 
proaches, the desired accuracy of a DMRG calculation 
can be established in advance if we take into account 
the dynamic change of the reduced density matrix of the 
subsystem. Within our approach, described in the next 
section, we will be able to show that if the number of 
block states is adjusted dynamically, a linear relationship 
obtains between the logarithm of the real error and the 
truncation error, which, in turn, can be used to derive a 
novel method to extrapolate to the full CI result. 

Our main goal in this paper is to determine the accu- 
racy of fc-DMRG in quantum chemistry and show that 
the algorithm converges to the accuracy that was set up 
in advance of the calculation. We have, therefore, car- 
ried out a detailed DMRG study of CH 2 , H 2 and F 2 
molecules with various number of orbitals each repre- 
senting different test cases. We have also addressed prob- 
lems related to the initial block-state configuration that 
arise within the framework of fc-DMRG. Since the focus 
of the paper is on the dynamic scaling of the density ma- 
trix and parameters of DMRG, we recall only those main 



definitions and formulas in this paper that are relevant 
to the question and not well known in quantum chem- 
istry. Therefore, details of our numerical procedure and 
developments will be published elsewhere. Although we 
have analyzed the general trend of the numerical error 
of fc-DMRG through quantum-chemical calculations, our 
results can be generally applied to other quantum system 
as well. 

The setup of the paper is as follows. In Sec. II we 
briefly describe the main steps of DMRG and recall the 
main sources of the numerical error. Sec. Ill is devoted 
to the details of the numerical procedure used to deter- 
mine the dynamic scaling behavior of the density matrix 
and to the problems that appear in the context of quan- 
tum chemistry. Sec. IV contains the numerical results 
and analysis of the observed trends of the numerical er- 
ror. The summary of our conclusions and a few general 
comments about the algorithm is presented in Sec. V. 

II. BACKGROUND OF THE NUMERICAL 
ERROR 

Detailed description of the_DMRG algorithm can be 
found in the original papersErfj and its application in the 
context of quantum chemistry is summarized in two re- 
cently published papersaD. Therefore, we present only 
the most important formulas and definitions that are rel- 
evant to the question of accuracy. 

The main purpose of DMRG is to treat the electron- 
electron correlation in a rigorous way which allows the 
minimization of the energy and calculation of measur- 
able quantities. Since DMRG is a variational procedure, 
it always provides an upper bound for all the calculated 
quantities. In the context of quantum chemistry a one 
dimensional chain that is studied by DMRG is built up 
from the molecular orbitals that were obtained, e.g., in 
a Hartree-Fock calculation. The electron-electron cor- 
relation is taken into account by an iterative procedure 
that minimizes the Rayleigh quotient corresponding to 
the Hamiltonian describing the electronic structure of the 
molecule, given by 

H X! / '-' r - r '" ■ H Vijklc\ a c\ el c krT ,c ll7 (1) 

ija ijklaa' 

and thus determines the full CI wavefunction. In Eq. (|l|) 
Tij denotes the matrix elements of the one-particle 
Hamiltonian comprising kinetic energy and the external 
electric field of the nuclei, and Vijki stands for the ma- 
trix elements of the electron repulsion operator. In or- 
der to show what are the key concepts and parameters 
of the numerical renormalization procedure and what are 
those drawbacks which hinder the analytical study of the 
method we have included a brief overview of the renor- 
malization group methods. 
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A. Block renormalization group method (BRG) 

In order to determine the eigenvalue spectrum of the 
Hamiltonian corresponding to an infinite long quantum 
chain (in the context of quantum chemistry this means 
infinitely many orbitals) built up from quantum sites rep- 
resented by q basis states, blocks were formed from each 
of two adjacent sites, and the Hamiltonian was deter- 
mined on the new configuration as is shown on Fig. [|. 
First the Hamiltonian of the model is diagonalized for 
two sites and then the q lowest energy states are selected 
out of the q 2 states whereby the so called block site will 
represent the two-sites problem in the subsequent itera- 
tion step. Operators defined on the selected q basis states 
are obtained from the original site operators according to 
a renormalization procedure given by the equation 



a. BRG 



(2) 



where operator O is constructed from the selected q 
eigenfunctions of the two-sites problem. In order to re- 
tain the original structure of the Hamiltonian operator, 
on-site (in the figure labeled by h) and inter-site (denoted 
by A) coupling constants are renormalized as well, shown 
as b! and A'. In the subsequent step the q 2 -dimensional 
Hamiltonian operator is diagonalized for two adjacent 
block sites, and again q states with lowest energy are 
selected out for the block site that will represent four 
sites in the following iteration step. Since the structure 
of the original Hamiltonian operator is retained and the 
number of coupling constants is unchanged, changes of 
the coupling constants (flow equations) can be studied 
analytically. When subsequent iteration of the renormal- 
ization steps leaves the coupling constants unchanged, 
the algorithm has reached a fix point which represents 
the infinite length (thermodynamic) limit of the model. 



B. Wilson's renormalization group method 

Besides a few analytically solvable models it turned 
out that the BRG method can be used only numeri- 
cally and its systematically increasing inaccuracy hin- 
dered the application of the method. In 1975 Wilson 
introduced another procedure for the numerical renor- 
malization methoclEjUj, in which a quantum chain with 
finite length L is built up systematically from quantum 
sites represented by q basis states by keeping the size of 
the Hilbert space fixed as is shown on Fig. |l|. 
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FIG. 1. Schematic plot of the spin couplings in the BRG, 
Wilson's and DMRG renormalization methods. 

The main idea of the method was again to solve the 
Hamiltonian of the model for two sites and selecting q' 
lowest energy states out of the q 2 states where q' was 
increased systematically up to a maximum value during 
the first few iteration steps based on the energy spectrum 
and kept constant afterwards. Operators were renormal- 
ized according to Eq. (|J). The key difference of Wilson's 
method compared to BRG is that he did not retain the 
original structure of the Hamiltonian operator, but he an- 
alyzed the scaling behavior of the energy as a function of 
the chain length. Systematical application of the renor- 
malization procedure introduces new terms and coupling 
constants. However, many of them become irrelevant for 
longer chains and the method also drives the system into 
the fixpoint. The major drawback of the method is that 
since the structure of the Hamiltonian changes with in- 
creasing chain lengths, flow equations can not be defined 
and the method can not be studied analytically. 



C. Density matrix renormalization group method 
(DMRG) 

In spite of the powerful properties of Wilson's proce- 
dure, the numerical error of the method grew systemat- 
ically with increasing chain length, which drawback has 
led to the fact that longer chains could not be studied 
numerically. Besides the truncation of the Hilbert space 
through the renormalization procedure, the numerical er- 
ror had another main source. When an additional un- 
renormalized site was added to the block site, the cou- 
pling was taken into account only between the block site 
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and this new site. In each iteration step the problem was, 
therefore, reduced to an isolated two-sites problem with 
open boundary conditions. These observations has led 
S. R. White to construct a larger auxiliary system (su- 
perblock configuration) which contains an environment 
in addition to the original block site problem to take 
care of the boundary effects in a more reliable way, as 
is shown on Fig. According to the figure the structure 
of the superblock configuration is defined as Bl * *Br 
where Bl represents the block site, • the new site under 
consideration, the additional *Br configuration the en- 
vironment and Ml and Mr denotes the number of block 
states, respectively. In order to minimize the error in- 
troduced in the representation of the block state in the 
truncation process, S. R. White has constructed the O 
matrix using the eigenfunctions of the reduced density 
matrix of the subsystem Bl*- It has been recognized in 
different contexttHl that the reduced subsystem density 
matrix describes the interactions of two subsystems in a 
particularly efficient way. Using these two key ingredi- 
ents, a DMRG iteration step first includes the diagonal- 
ization of the Hamiltonian constructed on the superblock 
configuration to obtain the target state. The target state 
is chosen from the eigenvalue spectrum of the Hamilto- 
nian that we want to calculate. It can also be a linear 
combination or even an incoherent superposition of more 
eigenstates as well. If \I) and |J) denote basis states 
for Bl* and *Br, respectively, then the target state is 
written as 

M L *q,M R *q 

^Target = £ M IJ )> ( 3 ) 
I,J 

where ipi t j is determined by diagonalization of the su- 
perblock Hamiltonian. After the target state is obtained, 
the reduced density matrix of the Bl* subsystem 

Pi, i' = ^ ipi,J^i',J (4) 
J 

is diagonalized and the M eigenstates with largest eigen- 
values (uj a ) are selected to build up the O matrix. The 
site operators are renormalized according to Eq. (Q) . The 
error of the truncation procedure in the DMRG method 
can be measured by means of the deviation of the total 
weight of the selected states from unity which is defined 
as 

Ml 

TRE = 1-Y j u a . (5) 

a = l 

The initial Bl and Br configuration contain one site 
per block each, thus the superblock Hamiltonian is de- 
termined on q 4 basis states restricted to the conserved 
quantum numbers like the total spin or the number of 
electrons. In each iteration step the size of the chain is 
increased by two sites until the desired chain length is 
reached as is shown on Fig. 111. This procedure is the so 



called infinite lattice algorithm. In order to average out 
long-wavelength fluctuations, the superblock configura- 
tion is asymmetrised by increasing the size of Bl and 
decreasing the size of Br until the left block contains 
L — 3 sites and the right block one site. The same proce- 
dure is then carried out in the reverse way and when the 
configuration is symmetric again, the first sweep of the so 
called finite lattice algorithm ends. This procedure can 
be repeated infinitely many times and is usually stopped 
when the energy does not change within two subsequent 
sweeps. There is again a major difference between BRG 
and DMRG which makes the analytical study of the scal- 
ing behavior of the latter method very complicated: In 
the DMRG method the number of selected block states 
(M) is larger then q and the original structure of the 
Hamiltonian is not retained, thus flow equations of the 
coupling constants can not be determined. 

According to the two key ingredients of the method, 
the numerical error of the DMRG algorithm has basi- 
cally two independent components which are the trun- 
cation error and the environmental error. The first one 
is generated during the renormalization step due to the 
truncation of Hilbert space, while the environmental er- 
ror appears because the chain is built up from blocks and 
the l ong range interactions are cut off. As it was shown in 
Ref. |l3j using the finite lattice method, the environmen- 
tal error can be averaged out and finally there remains 
a linear relationship on a log-log scale between the real 
error and truncation error. 

The truncation error, on the other hand, strongly de- 
pends on the shape of the eigenvalue spectrum of the 
reduced subsystem density matrix and on the number of 
block states kept for the subsequent iteration step. It 
has also long been known that the structure of the den- 
sity matrix depends on the criticality of the model. For 
systems with finite energy gap and coherence length the 
density matrix eigenvalue spectra decays exponentially, 
while for critical models with infinite coherence length it 
has a power-law tail. Besides these, in case of analytically 
solvable models the structure of the eigenvalue spectra of 
the density matrix determines the energy spectrum of the 
model as it was shown in Ref. |l4|. 

In addition to all the points discussed above, the de- 
cay of the eigenvalue spectrum also changes as the target 
state gets closer to the exact solution. It is, therefore, ev- 
ident that selecting out the M most probable states with 
highest eigenvalues will be an insufficient condition to 
control the accuracy of the DMRG method. Instead, one 
has to take care of the dynamic changes of the spectrum 
of the density matrix and keep the truncation error be- 
low a given threshold. Since the structure of the density 
matrix represents the whole system as well, it naturally 
arises that the number of block states should be selected 
out in a way that the truncation error satisfies an initial 
condition that was introduced in advance of the calcula- 
tion. 
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D. QC-DMRG method 

In the context of quantum chemistry, a one dimen- 
sional chain containing L molecular orbitals is generated 
by ordering the orbitals employed to build up the multi- 
particle states with increasing energy or by other rules, 
analogous to k points in fc-DMRGcl. These molecular or- 
bitals are calculated by standard numerical methods of 
quantum chemistry. 

It worth to note that the optimal ordering of the 
orbitals in the chain is still an open field of research. 
Note that the initial chain length of the QC-DMRG is 
L from the very beginning and the block operators for 
the left and right blocks are generated by a "warm up" 
proceduretj instead of the infinite lattice algorithm. The 
effect of the electron-electron correlation is taken into ac- 
count by the systematic sweeps in the framework of the 
finite lattice algorithm. Since the overall performance of 
the QC-DMRG method differs from the real-space ver- 
sion, it is also expected that new problems arise due to 
the inaccuracy of the starting wave function. These will 
be also investigated in detail in the next section. 

The most straightforward procedure to represent the 
unrenormalized site operators is to define them on spin- 
orbital basis states, in which case q is equal to two. The 
phase operator is then taken care of automatically by the 
standard definition of fcrmion creating and annihilating 
operators. On the other hand, if orbitals from, e.g., a re- 
stricted Hartree-Fock (RHF) calculations are employed, 
it is possible to define a super site built up from the 
ordered tensor product of spin-down and spin-up basis 
states, in which case q is 4 and the phase factor must 
be explicitly taken care of. This method offers consider- 
able efficiency gains because in this way the chain is only 
half the size compared to an unrestricted HF (UHF) type 
formulation, using spin orbitals for each site. Thus the 
number of multiplications using quadratic auxiliary oper- 
ators during the superblock Hamiltonian diagonalization 
procedures is roughly reduced by a factor of 4 compared 
to the spin-orbital formulation. In our implementation 
we have built up the chain from super sites. 

III. CONTROLLING THE ABSOLUTE ERROR 
OF DMRG 

A. Dynamic adjustment of the number of block 
states 

In order to control the accuracy of the DMRG proce- 
dure, the selection of the multi-particle states of the su- 
perblock Hamiltonian which are used for renormalization 
is obviously the decisive issue. Keeping all states featur- 
ing eigenvalues of the subsystem reduced density matrix 
larger than a fixed parameter which we called DM cut 
during the renormalization procedure, the truncation er- 
ror can be as small as DM cutl but it can be larger if 



the integrated contribution of the neglected states is still 
significant. To avoid such problem we propose to adjust 
DM cut dynamically, thus the number of selected states 
is increased as long as the integrated weight of neglected 
states is larger then a maximum value Ti?£' max , which 
can be fixed at the beginning of the calculation. This en- 
ables us to set up the desired accuracy of the DMRG al- 
gorithm at the beginning of the calculation. The number 
of states will be adjusted by this protocol in a dynamical 
fashion, depending on the structure of the density matrix 
spectrum. 

Since the truncation error is not immediately con- 
nected to the error in energy, one can control only the 
relative error in this way. In order to control the ab- 
solute error in energy, TRE max should be scaled by the 
Hartree-Fock energy or by the energy value calculated 
by the DMRG method which usually has the same or- 
der of magnitude as the exact value even after the first 
few iterations. We then expect the relative error of the 
energy to converge to this scaled threshold within a few 
sweeps of the DMRG procedure. 

From technical point of view, dynamic selection of 
block states has another important advantage. In the 
standard DMRG calculation the number of block states 
is fixed. Using our dynamical adjustment, the largest 
number of block states required to guarantee a given 
truncation error develops, however, only close to the sym- 
metric configuration during the sweep. For most of the 
remaining steps the threshold Ti?_E max is reached with 
considerably smaller number of block states, leading to 
substantial gains in efficiency in the renormalization step 
and the construction of the next superblock Hamiltonian, 
when dynamic block state selection is used. 

Within the framework of our procedure, it is also ev- 
ident why previously developed extrapolation methods 
based on functions of the number of block states failed 
to estimate the scaling behavior of the error in a rigorous 
way. The value of M is only one of the factors that de- 
termines the largest value of the truncation error during 
a full sweep. Using it exclusively, changes of the density 
matrix are not taken care of. Thus it is almost impossible 
to derive a reliable formula to estimate the real error as 
a function of the number of block states for the general 
case. 



B. Initial condition for the number of block states 

The straightforward application of dynamical control 
of DM cut during the first few sweeps is complicated by 
the fact that there is a major difference between the wave 
function of a given chain length generated by the infinite 
lattice algorithm of the real space version and that of gen- 
erated by ordering the orbitals in the case of A:-DMRG. In 
the first method, the wave function of the target state is 
always very close to the one which is obtained after sev- 
eral sweeps of the finite lattice method; however this is 
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not true in general for the momentum space version when 
the wave function strongly depends on the ordering of the 
orbitals. For example, it typically happens that during 
the first few steps the density matrix eigenvalue spectra 
will have very few states with large eigenvalues and many 
states with almost zero weight. In this case, the number 
of selected states will be cut drastically, which will limit 
seriously the size of the Hilbert space in the subsequent 
iterations, causing the algorithm being trapped in a local 
minimum. This situation happens in other optimization 
methods as well, and for example in the case of simulated 
annealing the so called adiabatic heating is used to move 
the algorithm out from the attractor of a local minimum. 
In the context of DMRG the introduction of virtual states 
is required in this situation, which means that we keep 
also those states that had almost zero eigenvalue up to 
a fixed number that we called Af m j n during the first two 
sweeps. Usually after the first sweep the decay of the 
density matrix spectrum becomes smooth and it changes 
dynamically as the target state gets closer to the exact 
one. 



C. New criteria for convergence and extrapolation of 
the FCI energy 

Up to now the condition for the number of sweeps 
was determined in an empirical way, using the condi- 
tion that the algorithm is stopped when the energy value 
obtained by two subsequent sweep does not change any 
more. Within the framework of DBSS we have a new 
criterion for the convergence. We have found that after 
convergence not only the energy value remains stable, 
but also the eigenvalue spectrum of the density matrix 
and thus the block states selected out by the algorithm 
for a given Bl • *Br configuration are the same during 
all subsequent sweeps. Although all subsequent sweeps 
leave the density matrix unchanged, still a fix point is 
not obtained since the structure of the density matrix 
and thus the truncation error and the relative error af- 
ter convergence can slightly change (but within the same 
order of magnitude) depending on the initial condition, 
for example on different M m j n . On the other hand, we 
can treat the energy values obtained for various TRE max 
values as points on a flow equation that converge to the 
fix point at the end^-which is the FCI energy. Based on 
our previous results£3 and those presented in the next 
section we can extrapolate to the FCI energy using the 
equation 

log ~ FCI = a * \og(TRE) + 6, (6) 
Efci 

where a, 6, Epci are parameters determined from the fit 
of the numerical result. As discussed below, our numeri- 
cal results show that the value of a is close to one. 



D. Error of the excited states due to the inaccuracy 
of starting block states 

Besides the problem of the initial structure of the den- 
sity matrix there is another difficulty which stems from 
the inaccuracy of the starting block wave functions. By 
contrast to the infinite lattice method when the tar- 
get state always remains in the same spin symmetry 
or changes sign periodically as a function of the chain 
lengthllj, the symmetry of the target state depends on 
the initial ordering in the case of /c-DMRG. This can 
lead to a major error, because the DMRG algorithm can 
lose the target state if its symmetry changes during the 
first few sweeps. It can happen that for example target- 
ing the second level the coefficients of the wave function 
of the ground state and excited states will mix and the 
spin symmetry of the target state changes randomly, and 
thus the energetically lowest level will be lost and the 
third level will become the target state. 

E. Introduction of local symmetry operators 

In order to avoid the random change of the spin sym- 
metry we have introduced partial spin adaption making 
sure that the permutational symmetry of the spins is 
odd for even S and even for odd S, which implements 
the spin reversal operator that flips the spins along the 
z-directions as it was shown in Ref. EM In case of k- 
DMRG the starting block wave function is constructed 
in a way that it contains basis states with N up and Ndown 
quantum numbers, thus fixing m s , and their symmetric 
components (i. e. , states with — m s ) as well. During the 
renormalization procedure a state and its partner belong 
to same eigenvalue of the density matrix, thus the dy- 
namic selection rule automatically ensures that both of 
them are kept. It worth to note, that this is not the 
full adaption of S 2 symmetry, which would be clearly be 
desirable, but more complicated to achieve in the frame- 
work of DMRG. Thus, components of the singlet and 
quintets levels can still mix, but it is not a problem since 
they are usually well separated. Application the corre- 
sponding spin reversal operator effectively ensures that 
the target state will remain in the spin symmetry sector 
that was fixed at the beginning of the calculation. 

From technical point of view, this has the additional 
advantage that one needs to target only the first level in 
both spin symmetry sectors, which always requires less 
block states to achieve a given accuracy. In addition, the 
number of auxiliary operators needed during the diago- 
nalization of the superblock Hamiltonian is decreased by 
a factor of two which doubles the speed of the algorithm. 
For the half-filled case the particle-hole symmetry oper- 
ator can be introduced in the same way. Details of the 
numerical procedure will be published elsewhere. 
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F. Error of the expectation value of one- and two 
particle operators 

The expectation value of the one- and two electron 
operators can be calculated from the one particle density 
matrix according to 

(A)=TR{pA). (7) 

where A is a L by L matrix of operator for a first-order 
property (e.g., dipole moment) in the same representa- 
tion as the original Tjj and Vijki were. So with A = Tij 
Eq. [?] provides the kinetic energy of the FCI wave func- 
tion. Once the target state was obtained, the one particle 
reduced density matrix can be formed for any con- 
figuration as 

Pij = (* Target | ^ c\ a C ja | ^Target ) (8) 

where i, j denote sites in the left block. The one-particle 
density matrix for the right block is determined in a sim- 
ilar way. If i is in the left block and j in the right block, 
then is constructed from the one- particle operators 
of the two blocks. The latter case was used to calcu- 
late two-point correlation functions in real space DMRG 
and was shown that the error of the one- and two-point 
correlation function is larger by one or two orders of mag- 
nitude compared to the error of the ground state energy. 
Since the dynamic block state selection rule controls the 
accuracy of the ground state, it also ensures the same 
scaling behavior of the correlation functions as well. Be- 
sides that, the fluctuation of the error shown in Rcf. [l3| 
because of the fluctuation of the truncation error within 
a full sweep caused by to the constant value of M also 
diminishes. 

The two-particle reduced density matrix can be ob- 
tained in a similar way 

Tijkl = {^Target] 22 c\ a c\ a , C ka > C k i a \^ Target) , (9) 
era' 

where the four-operator term is decomposed into four 
independent terms depending on the distribution of the 
i,j,k,l indices along the chain making use. of the usual 
partially contracted operators of fc-DMRQj. 

IV. NUMERICAL RESULTS 

In order to study the performance of fc-DMRG in quan- 
tum chemistry, we followed a route similar to the one ssic 
used to study the accuracy of the real-space DMRGE3. 
We performed calculations on molecules with different 
properties for which DMRG is expected to possess dif- 
ferent scaling behavior. Thus we have carried out a de- 
tailed DMRG study of the absolute error of the energy 
as a function of the number of orbitals and the fraction 
of filled orbitals on molecules CH2, H2O, and F2. The 



Hartree-Fock orbitals in a given basis of Gaussian or- 
bitals were calculated, and the T,j and Viju matrix ele- 
ments were transformed to the Hartree-Fock basis using 
the MOLPRO program packageEj, which was also [-used 
for the calculation of the benchmark full-CI energiestZHiS. 

We used various basis sets and geometries for the 
molecules which we selected for benchmark calculations. 
The geometries, references to the basis sets employed, 
and results obtained in SCF calculations as well as full 
CI energies are detailed in Table |. The models employed 
for the water molecule have also been used in White's 
studyQ. We include these cases here in order to enable a 
direct comparison with previous work. A more interest- 
ing test case was to study the CH2 molecule, for which 
we report energies for the triplet ground state as well 
as for the first excited (singlet) state. Hartree-Fock or- 
bitals of the closed-shell singlet configuration were em- 
ployed in all calculations on CH2. Calculations of the 
FCI energy of the triplet state were carried out in both 
the m s = and m s = ±1 spin sectors. In order to show 
that the relative error scales to TRE maK independently 
of the fraction of filled orbitals we have studied the half- 
filled chains by calculating the ground state of F2 with 
14 electrons and 14 orbitals (freezing the fluorine Is or- 
bitals and discarding the two highest virtual orbitals) and 
with 18 electrons and 18 orbitals. The latter calculation 
provides evidence that QC-DMRG is capable to provide 
cutting-edge CASSCF calculations with the potential to 
push their limits to active spaces well beyond a size which 
is feasible nowadays by standard methods. 

A. Dynamic selection of Block states 

QC-DMRG calculations on the water molecule demon- 
strate the dynamic selection of block states. In the first 
two panels of Fig. || we have plotted the number of block 
states which were selected in a calculation correlating 10 
electrons in 14 orbitals by means of QC-DMRG, starting 
with different values of Af m j n : 
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the region of the plateaus resulting in a very fast transver- 
sal of these regions. In addition, the maximum values of 
Ml and Mr occur at different iteration steps, thus for a 
given superblock configuration we find that even if one 
of them is very large, the other is usually much smaller. 
These two facts, finally, optimize the computational time 
and memory resources within a full sweep of the method. 

In order to show the dynamic change of the structure 
of the reduced subsystem density matrix, we have plot- 
ted in Fig. H the eigenvalues of the reduced subsystem 
density matrix obtained at the symmetric configuration 
(left and right blocks contained 6 orbitals) from a calcu- 
lation of the F2 molecule represented by 14 electrons and 
14 orbitals. 



FIG. 2. The figure shows the dynamically selected num- 
ber of left- and right block states Bl, Br, respectively, for 
two values of the minimum threshold value M m i n = 16, 64 
and the relative error as a function of iteration obtained with 
A/ m in = 16,64, 164. In all cases the 10 electrons of the H2O 
molecule were correlated in the double-zeta water model with 
14 orbitals, and T7?i5 max = 10 -10 was set in advance of the 
calculations. 

The number of block states for the left and right 
blocks is denoted by Ml and Mr, respectively. In 
the third panel of Fig. @ we give the relative error 
{(Edmrg — Epci) I Epci) of the calculation as a func- 
tion of M m j n and the the iteration step. The value of 
TRE max was set to 10~ 10 in advance of the calculations. 
It is evident from the figure that the maximum number 
of block state does not depend on the prescribed mini- 
mum value (M m j n ), although it is reached faster for larger 
Mmin- In order to show that the converged value of the 
accuracy does not depend on the threshold value (once a 
large enough value was taken) we have also included the 
result obtained with M m ; n = 164. It can be seen in the 
figure that the relative error converges to TRE max in all 
cases, but the speed of convergence strongly depends on 
M m in- In order to show that the QC-DMRG algorithm is 
trapped in a local minimum if M m i n is chosen too small, 
we carried out calculations with Af m ; n = 4, 8 and found 
indeed the number of block states being hindered to grow 
up. Similar test calculations on longer chains indicated 
that a larger value of M m i n = 64-100 is needed, thus we 
suggest that in order to to avoid problems related to lo- 
cal attractors and to obtain a faster performance a value 
of M m j n no less than 150-200 should be taken for longer 
chains. 

Investigating the scaling of the relative error shown in 
the third panel of Fig. one can find long plateaus where 
the accuracy of the method is not improved. In the usual 
DMRG calculations going through such plateaus costs al- 
most the same amount of time as calculating the region 
where the error drops significantly. By contrast, it can be 
seen on the figure that the minimum value of M occurs in 




FIG. 3. The figure shows the eigenvalue spectrum of the re- 
duced subsystem density matrix obtained for the F2 molecule 
after the end of the sweeps(S') of the finite lattice method. In 
the legend we have also included the number of selected block 
states (Ml, Mr) as a function of sweeps. 

There are several conclusions that one can draw from 
the figure. First of all, the density matrix spectrum de- 
cays very rapidly during the first few sweeps (S=0 is part 
of the "warm up" procedure ) which clearly implies the 
requirement of introduction of virtual states. On the 
other hand, as the target state gets closer to the FCI 
limit, the fraction of eigenvalues larger then 10 -15 in- 
creases significantly. It can be seen from the figure that 
the decay of the spectrum can be fitted by a linear line 
on a semilogarithmic scale for the largest eigenvalues, 
thus the density matrix spectrum decays exponentially, 
where the slope is related to the finite coherence length 
of the model. On the other hand, the slope of the line 
changes as a function of sweeps until the algorithm con- 
verges. Once the relative error converged to TRE max 
(which means for S > 7 in the case at hand) the slope of 
the decay remains the same, and this is the reason why 
the number of selected block states are the same for the 
subsequent sweeps. It worth to note that since the decay 
of the density matrix can be fitted by a straight line in 
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this model, the truncation error can be estimated as a 
function of the block states. However, in order to obtain 
a rigorous scaling behavior of the error as a function of 
block states one has to include the change of the slope as 
well, which in general strongly depends on the static and 
dynamic correlations of the models. 

B. Relationship between the relative error and 

In order to test that the relative error converges to a 
given value of TRE max we have ran independent calcu- 
lations for all the test molecules by adjusting TRE max 
from 10~ 3 up to 10~ n . The relative error of the first 
excited state obtained for the CH2 molecule with 6 elec- 
trons and 13 orbitals using M m ; n = 32 as a function of 
the iteration step and TRE max is shown on Figure [|. It 
can be seen on Fig. [|a. that the relative error of the first 
excited state also 
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FIG. 4. Calculations for CH2 with L = 13 sites shows the 
relationship between relative error and Ti?U max . The straight 
line is the result of the fit. GS denotes the triplet gound state, 
and 1XS denotes the first excited (singlet) ground state. 

converges to the values of TRE max set up in advance of 
the calculations. The converged value of the relative er- 
ror as function of TRE max for the first excited as well as 
for the ground state is plotted on Fig. |Jb. It is clear from 
the figure that there is a linear relationship between the 
converged value of the relative error and the truncation 
error, the fitted slope being 0.98. Fitting our results ob- 
tained for the various tests cases also with different M m j n 
values, we have found that the slope was always between 
0.95 and 1.1. Calculations in the m s = ±1 spin sec- 
tors provided a faster convergence for the ground state, 
as expected. The residual splitting of the m s = 0, 1, —1 
components of the triplet level was as low as 10~ 12 a.u. 

Calculations performed on the other test molecules 
with different number of basis states and for various 



values of M m - ln showed that the relative error scales to 
TRE max independently of the number of orbitals, frac- 
tion of filled orbitals and the threshold level of the num- 
ber of block states. Of course, the convergence gets 
slower for longer chain lengths and we usually needed 
6-8 sweeps to gain an absolute accuracy of 10~ 4 a.u. in 
the case of the CH2 molecule calculated with 57 orbitals. 

From the technical point of view, one can start a 
DMRG calculation by setting TRE max to 10 -3 and when 
the algorithm has converged (the energy is unchanged, 
the number of block states are unchanged, the slope of 
the density matrix remains the same ) TRE max can be 
adjusted by an order of magnitude until the desired max- 
imum value of the accuracy is reached. Using the calcu- 
lated energy values and the truncation error obtained 
for various values of TRE max (which is slightly below 
TRE max ) the FCI energy can be estimated by Eq. (||). 
This equation contains three free parameters {Epci 1 a, b) 
to determine from the fit. However, based on our results 
we can set the parameter a to one. We have found that 
one can gain one to three orders of magnitude improve- 
ments in the error of the correlation energy by the ex- 
trapolation method and that fixing the parameter a to 
one always provides an upper bound. In order to obtain 
a more accurate fit one needs more data points, thus 
TRE max should be adjusted in even smaller steps, es- 
pecially, if the calculations are carried out only up to 
a relative accuracy of 10~ 5 , but we have not done such 
analysis yet. 

In case of solid state physics, chains with various 
lengths are calculated and the thermodynamic limit is ex- 
trapolated by the the so-called finite-size scaling method. 
Using our procedure one can improve the energy values 
obtained for a given length L by two to three orders 
of magnitude, thus the overall performance of finite-size 
scaling procedure can be improved significantly. 

C. Scaling of the number of block states 

As we have shown, the number of block states depends 
on the structure of the reduced density matrix spectrum. 
Thus it is not possible to determine the scaling behavior 
of the maximum number of block states as a function of 
the number of orbitals and the fraction of filled orbitals 
in a rigorous way. On the other hand, in order to present 
a rough indication of computational resources used dur- 
ing our calculations we have we have collected the values 
of the maximum number of block states selected dynam- 
ically by our the method in Table ||. 

D. Other factors that affects the accuracy 

It is important to note that our scaling results are ob- 
tained only for a proper ordering of the orbitals in the 
initial chain. We have found that for some cases the 
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accuracy can be improved significantly if the HF levels 
were ordered with increasing energy (labeled by Ordi 
on Fig.||) while for other cases we had to "mirror" the 
chains and placed orbitals occupied in the HF configu- 
ration to the center of the chain (labeled by Ordi). A 
non-optimal ordering can in fact lead the method to be 
trapped by a local minimum. This situation is shown 
explicitly in Fig. ^| indicated by Ord\ . Even if M m i n was 
almost tripled, the relative 
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FIG. 5. The figure shows that an incorrect ordering can 
drive DMRG to a local minimum. 

error converged to the same local minimum, which, on 
the other hand, also supports our previous statements 
that M m ; n does not affect the final convergence signifi- 
cantly Changing the ordering, we have found for Ordi 
that the algorithm has always converged to the value of 
TRE max . Studying the optimaLprdering can be a ma- 
jor field of research. Chan et alu has already suggested 
a procedure to optimize the ordering in a recently pub- 
lished paper. We have not analyzed their solution yet. 

Another field of study can be the optimization of the 
superblock configuration. In our largest calculations for 
the half-filled case (18 electron in 18 orbitals) the number 
of selected block states grew up to 1500-1800 with sizes 
of the Hilbert space of the superblock configuration in- 
creasing beyond 1.000.000. In order to decrease the size 
of the Hilbert space we have derived an alternate proto- 
col, modifying the superblock configuration as Bl • Br 
in a similar way as was done by Xiang. We found a con- 
siderably worse performance for these calculations. We 
believe that in the future it can worthwhile analyzing 
the speed of convergence for various superblock configu- 
rations. 



V. SUMMARY 

We have applied the momentum-space version of 
DMRG in quantum chemistry in order to study the accu- 
racy of the method. Analyzing the eigenvalue spectrum 
of the reduced density matrix and based on our previous 
results obtained for real space DMRG, we have shown 
that it is possible to set up the accuracy of the method in 
advance of the calculation by dynamically controlling the 
truncation error and the number of block states. We have 
carried out a detailed QC-DMRG study of the molecules 
H2O, CH2, and F2 obtained with various basis sets in 
order to show that the relative error scales to the max- 
imum threshold value of the truncation error that was 
fixed in advance of the calculation. We found that the 
linear relationship between the logarithm of the relative 
error and the logarithm of the maximum value of the 
truncation error is independent of the number of orbitals 
and the fraction of filled orbitals for the cases consid- 
ered. Based on these results we have presented a novel 
approach to extrapolate the FCI energy, a method which 
could also improve the accuracy of the finite-size scaling 
method when fc-DMRG is applied in solid state physics. 
We have addressed new problems related to the inaccu- 
racy of starting block state configuration and presented 
solutions to achieve faster convergence and better stabil- 
ity of the target state. 

The maximum number of block states that the algo- 
rithm selected out in the dynamic fashion was in the 
range of 1500-2000, the largest size of the Hilbert space 
related to the superblock Hamiltonian was 800.000- 
1.200.000, and the longest chains that we have studied 
contained 57 sites. 

Although momentum is not a good quantum number 
if fc-DMRG is applied in quantum chemistry, there are 
still a few remarks which might indicate why QC-DMRG 
method can work well in the field: 

• In most of the cases the calculations are carried out 
in the small U limit known to converge fast. 

• The number of electrons is fixed for a given molecule. 
Therefore, doubling the length of the system will not im- 
ply in general keeping the fraction of filled orbitals fixed. 
Thus calculating a molecule with more basis states would 
mean longer chains but with lower filling value which usu- 
ally has a better convergence. 

• The practical use of DMRG in quantum chemistry 
can open a route to active spaces well beyond today's 
limits, yielding complete active space self-consistent field 
(CASSCF) solutions with a relative error of the correla- 
tion energy of the order of 10~ 4 to 10~ 5 . This can be 
realized by a few thousand block states, which also ex- 
pected to hold for longer chains as well. Therefore, we 
believe that 3000-4000 block state will provide satisfac- 
tory results for all the chain lengths and fillings which 
are of interest in the immediate future. 

• Although the structure of the Hamiltonian is very 
complicated, it is decomposed into several parts. This 
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means that during the diagonalization step each compo- 
nent of the Hamiltonian can be applied on the wavefunc- 
tion independently, therefore, the method is an excellent 
candidate for parallel computers. 

Our source code was written in the the framework of 
the Matlab programming environment and the C++ code 
as well as the standalone code was produced by the Mat- 
lab compiler. Most of our numerical calculations were 
carried out on Athlon XP 1800+ processors under Linux 
and in some cases on a SGI 3000 machine of the local 
computer center. For the largest calculations comprising 
M =1700-2000 block states (F 2 18/18) the program re- 
quired 200-500 MB of RAM, running about 50-60 hours 
on Athlon XP 1800+ processor to achieve 10 -4 -10 -5 a.u. 
absolute accuracy. The scaling of computational time 
with the number of orbitals still can not be determined 
because of the development stage of our code, but as a 
rough indication it took some 15 hours for 8/24 chain and 
more than a week for the 6/57 case. The present stage 
of our code limited the number of block states around 
2000, however, solving a few technical points we expect 
that the feasible M can be increased significantly in the 
future. 
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reference 




split valencea 



(a.u.) 

1.84345 
1.88973 
2.02230 

2.02230 

2.02230 
2.64373 
2.68797 



(degrees) 

110.565 
104.500 
129.4667 

129.4667 

129.4667 



10 

8 
6 
6 
6 
6 
6 

14 

18 



14 
25 (24) 
14 (13) 
14 (13) 
24 (23) 
24 (23) 
58 (57) 
20 (14) 

18 



TABLE I. Geometries and benchmark energy values for 
the calculated molecules. The number of correlated orbitals 
is given in parentheses, unless it agrees with the total number 
of orbitals. 
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CH 2 


H 2 


F 2 


CH 2 


H 2 


F 2 


CH 2 


L 


6/13 


10/14 


14/14 


6/23 


8/24 


18/18 


6/57 


Filling 


0.230 


0.357 


0.500 


0.130 


0.166 


0.500 


0.052 


^EAba 


M max 


M max 




M max 


M max 


M max 


M max 




25 


40 


280 


150 


130 


170 


300 


10~ 3 


40 


60 


350 


280 


320 


520 


480 


1(T 4 


100 


140 


800 


370 


440 


1100 


620 


10- 5 


160 


300 


1450 


580 


650 


1800 




10~ 6 


230 


420 




670 


820 






10- 7 


300 


530 




720 








10~ 8 


360 


650 




880 








1(T 9 


420 


780 













TABLE II. The maximum number of the block states se- 
lected out dynamically by DMRG to reach a given value of 
absolute accuracy. The second row contains the number of 
electrons and orbitals of each test calculations and below the 
fraction of filled orbitals is listed. 
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